Relationships among seismic, electrical and reservoir properties are often exploited in geophysical exploration to model geophysical properties of subsurface regions, e.g., where data from seismic and/or electromagnetic surveys are used to predict a range of features of a subsurface region. The predicted geophysical features are then used for various exploration decisions, e.g., the number of wells to drill, the type of well(s) to drill, and optimal well location to recover resource(s) from a reservoir.
Seismic properties of a subsurface region are those properties that directly determine the reflection and transmission of seismic waves by the subsurface, and together define at least the compressional wave velocity, shear wave velocity, and density of the subsurface region. It is often more convenient to express the seismic properties of a subsurface in terms of elastic properties, such as bulk modulus and shear modulus (also called the elastic moduli). Various functions of the velocities and density of the subsurface can also be equivalently used to express seismic properties, including: bulk modulus, Poisson's ratio, Vp/Vs ratio, P-wave modulus, impedance, and Lamé parameters. Seismic properties may also include, for example, anisotropy and attenuation. Seismic wave velocities may also vary with the frequency of the seismic wave, a phenomenon called dispersion.
Among the seismic properties, impedance is the product of seismic velocity and the density. Impedance, also called acoustic impedance and often symbolized by IP, will typically vary among different rock layers, e.g., opposing sides of an interface will have different impedances. The reflection coefficient of an interface generally depends on the contrast in acoustic impedance of the rock on either side of the interface. Specifically, the difference in acoustic impedance between rock layers affects the reflection coefficient. One geophysical modeling process for determining the impedance structure of a subsurface region based on recorded seismic reflection data is seismic inversion.
Seismic inversion techniques rely upon seismic reflection data, typically obtained through a seismic survey and analysis of the seismic data from the survey. Seismic reflection techniques are typically based on the generation of seismic waves in the earth's surface, through the use of one or more seismic sources, e.g., dynamite, air guns, vibrators, and the recording and analysis of the portions of these waves that get reflected at the boundaries between the earth's layers. FIGS. 1A-1B are views of convolutional models for seismograms generated from primary reflections at one or more boundaries between two or more media. Referring to FIG. 1A, a single boundary model 100 shows that at a given boundary between two media, the amplitude (strength) of the reflected wave is proportional to the amplitude of the incident wave and a quantity called a reflection coefficient. The value of the reflection coefficient depends on the elastic parameters of the two media, and for normal incidence it is given by equation (1). The seismic trace for this case contains a single pulse, whose shape is that of the seismic wavelet.
The reflection coefficient, for normal incidence (rays perpendicular to the reflecting interface), is defined as:R=(IP2−IP1)/(IP2+IP1)  (1)In equation (1), R is the reflection coefficient and the quantities IP1 and IP2 are called compressional impedances.
The terms P-impedance and acoustic impedance are also commonly used to describe the same quantities. For example, compressional impedance is defined as the product of density and compressional (P-wave) velocity:IP=ρVP  (2)In this equation ρ is density and VP is the P-wave velocity. In equation (1), IP1 and IP2 are the compressional impedances of the layers above and below the reflective interface, respectively. For a large number of reflecting boundaries, the recorded seismic reflection response is the sum of the responses for the different boundaries.
Referring to FIG. 1B, the multiple boundary model 150 shows that a reflection event is typically recorded on every seismic trace at any given time. The recorded seismogram for the multiple boundary reflection configuration can then be thought of as a reflectivity time series, e.g., that is symbolized by r(t) and based on an impedance profile IP(t). If multiple reflections are ignored, and the pulse generated by the seismic acquisition system is a simple spike, the recorded seismic trace is composed of a sequence of reflectivity spikes, with the size of each of them calculated based on equations (1) and (2).
However, the incident seismic wave is typically not a simple spike, but a broader waveform, called the seismic wavelet w(t). In this case, the recorded seismogram is not be r(t). Instead, every spike is replaced by an appropriately scaled version of the seismic wavelet, and the results added. When the reflecting medium contains multiple reflecting boundaries, the resulting seismic trace is further evaluated by calculating the convolution of the seismic wavelet and the reflectivity time series. The reflectivity time series is a sequence of spikes, each of them generated by a single boundary, according to equation (1). The mathematical operation that combines the reflectivity time series r(t) and the wavelet w(t) in the manner just described is convolution:s(t)=r(t)*w(t)  (3)where the symbol * denotes the operation of convolution in equation (3). In equation (3), the recorded seismogram s(t) is calculated as the convolution of the reflectivity series r(t) and the wavelet w(t). Equation (3) expresses what is typically referred to as the convolutional model of reflection seismology.
Assuming continuous recording of seismic reflections, the equation for calculating the normal-incidence reflection coefficient (equation (1)), can be generalized to the following expression:r(t)=(dIP(t)/dt)/(2IP(t))  (4)
In equation (4), IP(t) represents the impedance value for a layer at a depth such that the reflection from the layer is recorded at a time t. The operator d/dt represents the derivative with respect to time. An exemplary seismic inversion problem from normal-incidence seismic data amounts to solving equations (3) and (4) to determine the impedance function IP(t), and assuming knowledge of the recorded seismic data s(t) and the seismic wavelet w(t). In the limit when the time interval between recorded spikes is very small, one can consider the reflectivity series as a continuous function of time, whose relationship to impedance, for normal incidence, is given by equation (4). For non-normal incidence the calculation of the reflection coefficients is modified, but the convolutional model, as described here for primary reflections only, remains valid.
Estimation of the seismic wavelet w(t) can be achieved by making use of well log data. When a well is available and appropriate sonic and density well logs have been collected, the impedance IP(t) and reflectivity r(t) are known. Equation (3) can then be used to solve for w(t), given r(t) and the seismic trace s(t). For this estimation to work adequately, an accurate correlation usually needs to be established between subsurface information at the well and the seismic events. The term “well tie” is commonly used to describe the process of establishing this correlation. Accordingly, accurate well ties are a prerequisite for most inversion methods.
The aforementioned concepts can also be generalized to the case where the recorded reflections correspond to larger angles between the incident and reflected wave propagation paths, e.g., oblique or non-normal incidence case. For such situations the convolutional model equation (3) is still valid, but the expression for the reflection coefficient equation (4) is replaced with a more complicated expression, e.g., containing additional elastic parameters, such as shear-wave velocity.
Various seismic inversion techniques based on the convolutional models have been applied in common practice. Two recently developed seismic inversion techniques that are implemented as simple modifications of the frequency spectrum are Coloured inversion and Spectral Shaping inversion. These seismic inversion techniques are further described in Lancaster, S., and Whitcombe, D., 2000, “Fast Track “Coloured” Inversion,” Expanded Abstracts, 70th SEG Annual Meeting, Calgary, 1572-1575; and Lazaratos, S., 2006, “Spectral Shaping Inversion For Elastic And Rock Property Estimation,” Research Disclosure, Issue 511, November 2006.
Referring to FIG. 2, while the two techniques differ in their implementation, both inversion techniques are similar conceptually. For example, impedance estimation is performed by the combination of a phase rotation (−90°) and a spectral shaping operation applied to seismic data. Prior to the application of the phase rotation and the spectral shaping operation, the seismic data is typically converted to zero-phase, e.g., for zero-phase data all frequency components of the seismic wavelet are synchronized and combined to produce a wavelet that is symmetric around the wavelet peak. Coloured inversion assumes a log amplitude spectra follows an exponential law, while spectral shaping inversion (Lazaratos) does not require this assumption. In addition, coloured inversion is strictly a zero-offset inversion. Spectral shaping inversion also provides added benefits of being useful in generating estimates of both elastic and rock properties.
The spectral shaping operation is implemented by the application of a filter that reshapes the original seismic spectrum to make the seismic spectrum similar to the average spectrum of well logs recorded at wells in the subsurface region. Referring to FIG. 2, a graphical view 200 demonstrates how the spectral shaping filters significantly amplify the energy in the low-frequency part of the seismic spectrum. Average local well log 220 and original seismic frequency 240 spectra are significantly different even over the range of frequencies for which the signal-to-noise ratio of the data is positive. Spectral shaping reshapes the original spectrum to make it similar to the log spectrum. The resulting frequency spectrum is the shaped seismic spectrum 260. The shaping operation implies significant amplification of the low-frequency energy, as seen in FIG. 2.
Lazaratos [2006] provides a mathematical derivation demonstrating that, under assumptions that are commonly satisfied, the spectral shaping procedure highlighted above provides an estimate of the impedance, solving equations (3) and (4). For example, based on the convolutional model established above, a seismic trace can be expressed by the convolution equation (5):
                              s          ⁡                      (            t            )                          =                                            w              ⁡                              (                t                )                                      *                          r              ⁡                              (                t                )                                              =                                                    w                ⁡                                  (                  t                  )                                            *                                                Δ                  ⁢                                                                          ⁢                                                            I                      P                                        ⁡                                          (                      t                      )                                                                                        2                  ⁢                                                            I                      P                                        ⁡                                          (                      t                      )                                                                                            =                                                            Δ                  ⁢                                                                          ⁢                  t                                2                            ⁢                              w                ⁡                                  (                  t                  )                                            *                              1                                                      I                    P                                    ⁡                                      (                    t                    )                                                              ⁢                                                ⅆ                                                            I                      P                                        ⁡                                          (                      t                      )                                                                                        ⅆ                  t                                                                                        (        5        )            
In the above expression, and hereinafter, the following notation convention is used to describe one or more of the following features:                s(t), S(f) seismic trace and its Fourier transform        Squad(f) Fourier transform of quadrature trace        w(t), W(f) wavelet and its Fourier transform        r(t) reflectivity        IP(t), IP(f) P-impedance and its Fourier transform         IP lowpass filtered P-impedance        Δt the sampling rate        
The term IP(t) in the denominator can be replaced by a very slowly changing function, which just contains the trend in IP. In practice, such a function can be generated by lowpass filtering IP, to maintain frequencies at the very low end of the spectrum (e.g. 0-2 Hz). This low-frequency term can then be treated as a simple multiplier and moved to the left of the convolution operator. The convolution equation then becomes (equation (6)):
                              s          ⁡                      (            t            )                          =                                                            Δ                ⁢                                                                  ⁢                t                            2                        ⁢                          w              ⁡                              (                t                )                                      *                          1                                                                    I                    P                                    _                                ⁡                                  (                  t                  )                                                      ⁢                                          ⅆ                                                      I                    P                                    ⁡                                      (                    t                    )                                                                              ⅆ                t                                              =                                                    Δ                ⁢                                                                  ⁢                t                                            2                ⁢                                                      I                    P                                    _                                                      ⁢                          w              ⁡                              (                t                )                                      *                                          ⅆ                                                      I                    P                                    ⁡                                      (                    t                    )                                                                              ⅆ                t                                                                        (        6        )            
A weak-scattering assumption, stated as follows, is relied upon to mathematically show the ability to transform the convolution equation from its original form to the one given in equation (6). P-impedance can be decomposed into a slowly changing background part, e.g., low-frequency trend, frequencies well below the seismic bandwidth, and a higher-frequency perturbation part including changes in the seismic bandwidth and above. Accordingly, (i) the perturbation should be weak relative to the background, and (ii) the background is essentially constant within the length of the seismic wavelet. Based on numerous observations supporting these conclusions, transforming equation (6) to the frequency domain results in equation (7):
                                          S            quad                    ⁡                      (            f            )                          =                                            πΔ              ⁢                                                          ⁢              t                                                      I                P                            _                                ⁢                      fW            ⁡                          (              f              )                                ⁢                                    I              P                        ⁡                          (              f              )                                                          (        7        )            
Averaging for several wells (using < > to signify the averaging operation), results in equation (8):
                              〈                                    S              quad                        ⁡                          (              f              )                                〉                =                                            πΔ              ⁢                                                          ⁢              t                                                      I                P                            _                                ⁢                      fW            ⁡                          (              f              )                                ⁢                      〈                                          I                P                            ⁡                              (                f                )                                      〉                                              (        8        )            where it is assumed that the seismic wavelet is constant for the area within which the wells are located.
By definition, the shaping filter's frequency response is the ratio of the average log spectrum and the average seismic spectrum, as seen in equation (9):
                              Shaping          ⁢                                          ⁢          Filter          ⁢                      :                    ⁢                                          ⁢                      H            ⁡                          (              f              )                                      =                                            〈                                                I                  P                                ⁡                                  (                  f                  )                                            〉                                      〈                                                S                  quad                                ⁡                                  (                  f                  )                                            〉                                =                                                                      I                  P                                _                                            πΔ                ⁢                                                                  ⁢                t                                      ⁢                          1                              fW                ⁡                                  (                  f                  )                                                                                        (        9        )            and applying this to the seismic data results in equation (10):Shaped Seismic=H(f)Squad(f)=IP(f)  (10)
Seismic migration of seismic data is a correction technique involving rearrangement of seismic events, so that reflections are plotted at a true representation of their subsurface locations. Referring to FIG. 3, a graphical model 300 shows, on the original recorded data, reflections from dipping interfaces are recorded at surface positions that are not directly above the subsurface locations where the reflections take place. In addition, isolated point-like discontinuities in the subsurface (point scatterers) generate seismic events (diffractions) recorded over a large range of receivers, that can make the interpretation of seismic data confusing. Seismic velocity variations are one more reason the original recorded data provide only a distorted view of the subsurface geology. The seismic migration technique addresses the above issues and is therefore utilized in many seismic data processing sequences to accurately depict the structures and geometric configurations observed in seismic recordings as an analog of the geologic layers that gave rise to the seismic reflections.
The need to correctly position dipping reflectors is best seen in FIG. 3. The reflection pulse from point A generated from a source at S1 and recorded at a receiver also at S1 is plotted on the trace under S1, at point A′, which is selected such that the lengths of S1A and S1A′ are equal (assuming a constant-velocity subsurface for simplicity). Similarly, the reflection pulse from point B is plotted on the trace under S2, at point B′. The reflector segment AB is plotted at the erroneous lateral position A′B′ and has a dip smaller than AB's true dip. Migration is the correction technique that corrects such distortions. Before migration, the structures and geometric configurations observed in seismic recordings are typically not an accurate description of the geologic layers that gave rise to the seismic reflections.
Seismic inversion has traditionally been limited to applications where seismic inversion has been applied after migration as accurate well ties are typically required to estimate the seismic wavelet. Since the original “un-migrated” data forms an inaccurate structural image of the sub-surface, accurate well ties are typically established after migration. The present inventors have determined that there is a need for a seismic inversion technique that can be applied at various stages in a modeling process while still being computationally efficient and accurate when used in conjunction with a migration correction technique to model impedance of a subsurface region.